#include "mloop.h"

/* data structure for topology */
struct bond {
	struct atomname atom1, atom2;
	int free;
} *bond;
int M;

/* data structure for connect table */
struct connect {
	int n, list[4], free[4];
} connect[LOOPNUM];

void read_top()
{
	FILE *fp;
	char buf[256];
	int i;

	fp = fopen ("TOPOLOGY", "r");
	if (fp == NULL)
	{
		printf ("open TOPOLOGY failed\n");
		exit (-1);
	}
	M = 0;
	while (fgets (buf, 256, fp))
	{
		if (!strncmp (buf, "BOND", 4)) M++;
	}
	fclose (fp);
	bond = (struct bond *) malloc (M * sizeof (struct bond));
	fp = fopen ("TOPOLOGY", "r");
	i = 0;
	while (fgets (buf, 256, fp))
	{
		char resname[4];
		if (!strncmp (buf, "RESI", 4))
		{
			strncpy (resname, buf + 5, 3);
		}
		else if (!strncmp (buf, "BOND", 4))
		{
			strncpy (bond[i].atom1.atomname, buf + 5, 3);
			strncpy (bond[i].atom2.atomname, buf + 9, 3);
			strncpy (bond[i].atom1.resname, resname, 3);
			strncpy (bond[i].atom2.resname, resname, 3);
			sscanf (buf + 13, "%d", &bond[i].free);
			i++;
		}
	}
	fclose (fp);
}

void check_connect()
{
	int i, j, k, n;
	struct loopatom *atom;

	n = loop.loopnum;
	atom = loop.loopatom;
	for (i = 0; i < n; i++)
		connect[i].n = 0;

	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			if (i == j)
				continue;
			if (atom[i].resno - atom[j].resno > 1)
				continue;
			/*
			if (atom[j].resno - atom[i].resno > 1)
				break;
			*/
			if (atom[i].resno == atom[j].resno)
			{
				for (k = 0; k < M; k++)
				{
					if (strncmp (atom[i].name.resname, bond[k].atom1.resname, 3) == 0 && strncmp (atom[i].name.atomname, bond[k].atom1.atomname, 3) ==0 && strncmp (atom[j].name.atomname, bond[k].atom2.atomname, 3) == 0)
					{
						connect[i].list[connect[i].n] = j;
						connect[j].list[connect[j].n] = i;
						connect[i].free[connect[i].n] = bond[k].free;
						connect[j].free[connect[j].n] = bond[k].free;
						connect[i].n++;
						connect[j].n++;
						break;
					}
				}
			}
			else if (atom[j].resno - atom[i].resno == 1)
			{
				if ( strncmp (atom[i].name.atomname, "C  ", 3) == 0 && strncmp (atom[j].name.atomname, "N  ", 3) == 0 )
				{
					connect[i].list[connect[i].n] = j;
					connect[j].list[connect[j].n] = i;
					connect[i].free[connect[i].n] = 0;
					connect[j].free[connect[j].n] = 0;
					connect[i].n++;
					connect[j].n++;
				}
			}
		}
	}
}

void checktable(int i, int j, int k)
{
	int p, q, r;
	for (p = 0; p < connect[j].n; p++)
	{
		q = connect[j].list[p];
		r = connect[j].free[p];
		if (i == q || table[i][q])
			continue;
		if (k == 1)
		{
			table[i][q] = 1;
			checktable (i, q, k + 1);
		}
		else
		{
			table[i][q] = 1;
			if (r)
				continue;
			else
				checktable (i, q, k + 1);
		}
	}
}

void checktable14(int i)
{
	int j, k, l, p, q, r;

	for (p = 0; p < connect[i].n; p++)
	{
		j = connect[i].list[p];
		for (q = 0; q < connect[j].n; q++)
		{
			k = connect[j].list[q];
			for (r = 0; r < connect[k].n; r++)
			{
				l = connect[k].list[r];
				if (!table[i][l])
					table[i][l] = 2;
			}
		}
	}
}

void check_table()
{
	int i, j;

	for (i = 0; i < loop.loopnum; i++)
		for (j = 0; j < loop.loopnum; j++)
			table[i][j] = 0;

	for (i = 0; i < loop.loopnum; i++)
	{
		checktable (i, i, 1);
		checktable14 (i);
	}
}

void build_table()
{
	read_top ();
	check_connect ();
	check_table ();
}
